###################################################################################
# Replication files for: Negative References to Amicus Briefs in Judicial Reasoning
# Authors: Johan Lindholm, Daniel Naurin & Philipp Schroeder
# Published in: Journal of Law and Courts
# Contact: P.Schroeder@gsi.uni-muenchen.de

### File 1: This script provides the findings and figures presented in the main manuscript

library(arm)
library(rstanarm)
library(broom)
library(tidyverse)
library(stringr)
library(pscl)
library(texreg)
library(sandwich)
library(lmtest)
library(ggplot2)
library(boot)
library(parallel)
library(margins)
library(lavaan)
library(lavaan.survey)
library(semTable)
library(Zelig)
library(nnet)

rm(list=ls())
homeFolder <- Sys.getenv("HOME")


# read data ---------------------------------------------------------------

load("data/mentions_issue_level.RData")
load("data/mentions_submission_level.RData")
load("data/post_nice.RData")



# path analysis for figure 2 ----------------------------------------------

path1 <- '
  cjeu_pos_auto.factor ~ com_pos_auto.numerical + ag_pos_auto.numerical + ag_com_disagreement + z.ms_pos_auto0 + z.ms_pos_auto1
  mentioned_negative_dummy ~ cjeu_pos_auto.factor + z.ms_pos_auto0 + z.ms_pos_auto1 + ag_cjeu_disagreement + ag_com_disagreement +
    z.sources_case_law.ag + z.sources_primary.ag + z.sources_secondary.ag +
    z.duration_since_ag + panel_size +
    number_of_observations + legact_primary + legact_secondary +
    administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
    labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport +
    z.subject_matter.count + year_proceeding.factor
'

path1.fit <- sem(path1, data = mentions_issue_level)
summary(path3.fit)
round(fitMeasures(path1.fit, c("chisq", "cfi", "rmsea", "pvalue")), 2)

path2 <- '
  cjeu_pos_auto.factor ~ com_pos_auto.numerical + ag_pos_auto.numerical + ag_com_disagreement + z.net_ms_autonomy
  mentioned_negative_dummy ~ cjeu_pos_auto.factor + z.net_ms_autonomy + ag_cjeu_disagreement + ag_com_disagreement +
    z.sources_case_law.ag + z.sources_primary.ag + z.sources_secondary.ag +
    z.duration_since_ag + panel_size +
    number_of_observations + legact_primary + legact_secondary +
    administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
    labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport +
    z.subject_matter.count + year_proceeding.factor
'
path2.fit <- sem(path2, data = mentions_issue_level)
summary(path2.fit)

# fit statistics for table 6 in the manuscript's supplementary materials
round(fitMeasures(path2.fit, c("chisq", "cfi", "rmsea", "pvalue")), 2)

GAMMA1_coefs <- parameterEstimates(path1.fit, ci = TRUE, level = 0.95)[1:5,c(3:4,8:9)]
BETA1_coefs <- parameterEstimates(path1.fit, ci = TRUE, level = 0.95)[6:18,c(3:4,8:9)]

GAMMA2_coefs <- parameterEstimates(path2.fit, ci = TRUE, level = 0.95)[1:4,c(3:4,8:9)]
BETA2_coefs <- parameterEstimates(path2.fit, ci = TRUE, level = 0.95)[5:16,c(3:4,8:9)]


# plot coefficients for GAMMA and BETA matrices for figure 2
pdf(file = "figures/path_coefficients.pdf", height = 11, width = 10)
par(mar=c(5,18,3,3), mfrow = c(2,1))
# GAMMA
plot("", axes = FALSE,
     xlim = c(-0.5,1.5),
     ylim = c(0.75,6.25),
     xlab = "Coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Path analysis: CJEU position (Models 1 and 2)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-0.5,1.5,.25), cex.axis = 1.2)
axis(2,at = seq(6,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = c("Commission position",
               "AG position",
               "Conflicting positions: AG/Commission",
               "MS favouring no restrictions",
               "MS favouring restrictions",
               "Net position of MS"
     ))
abline(h = seq(6,1,-1), lty = 2, lwd = .5, col = "grey")
# path 1
points(GAMMA1_coefs[,2], seq(6,2,-1) + .1, pch = 19)
segments(GAMMA1_coefs[,3], seq(6,2,-1) + .1,
         GAMMA1_coefs[,4], seq(6,2,-1) + .1, lwd = 2)
# path 2
segments(GAMMA2_coefs[,3], c(seq(6,4,-1),1) -.1,
         GAMMA2_coefs[,4], c(seq(6,4,-1),1) - .1, lwd = 2)
points(GAMMA2_coefs[,2], c(seq(6,4,-1),1) - .1, pch = 24, bg = "white")
abline(v=0, lty = 2)
# BETA
plot("", axes = FALSE,
     xlim = c(-0.1,0.1),
     ylim = c(0.75,14.25),
     xlab = "Coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Path analysis: Negative reference (Models 1 and 2)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-0.1,0.1,.05), cex.axis = 1.2)
axis(2,at = seq(14,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = c("CJEU position",
               "MS favouring no restrictions",
               "MS favouring restrictions",
               "Net position of MS",
               "Conflicting positions: AG/CJEU",
               "Conflicting positions: AG/Commission",
               "Sources cited by AG: Case law",
               "Sources cited by AG: Primary law",
               "Sources cited by AG: Secondary law",
               "Deliberation time",
               "Panel size",
               "Number of observations",
               "Issue concerns primary EU law",
               "Issue concerns secondary EU law"))
abline(h = seq(14,1,-1), lty = 2, lwd = .5, col = "grey")
# path 1
points(BETA1_coefs[,2], c(seq(14,12,-1), seq(10,1,-1)) + .15, pch = 19)
segments(BETA1_coefs[,3], c(seq(14,12,-1), seq(10,1,-1)) + .15,
         BETA1_coefs[,4], c(seq(14,12,-1), seq(10,1,-1)) + .15, lwd = 2)
# path 2
segments(BETA2_coefs[,3], c(14, seq(11,1,-1)) - .15,
         BETA2_coefs[,4], c(14, seq(11,1,-1)) - .15, lwd = 2)
points(BETA2_coefs[,2], c(14, seq(11,1,-1)) - .15, pch = 24, bg = "white")
abline(v=0, lty = 2)
dev.off()



# logistic regression for figure 3 ----------------------------------------

# subset post nice treaty
m3 <- glm(mentioned_negative_dummy ~ ag_opinion_submitted + panel_size +
            number_of_observations + com_pos_auto.factor +
            z.ms_pos_auto0 + z.ms_pos_auto1 + 
            year_proceeding.factor + z.subject_matter.count + 
            legact_primary + legact_secondary +
            administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
            labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport,
          data = post_nice, family = binomial(link = "logit"))

m3coeffs_cl <- cluster_se(m3)
vcov_mat3 <- vcovCL(m3, cluster = ~iuropa_proceeding_id)

# check for multicollinearity
car::vif(m3)

# coefficient estimates
get_point_estimates <- function(model, se, var_index){
  bin <- matrix(nrow = 3, ncol = length(names(coef(model))))
  bin[1,] <- coef(model)
  bin[2,] <- coef(model) - 1.96 * se[,2]
  bin[3,] <- coef(model) + 1.96 * se[,2]
  bin <- as.data.frame(bin)
  names(bin) <- names(coef(model))
  
  # select coefficients
  out <- as.data.frame(matrix(nrow = 3, ncol = length(var_index)))
  names(out) <- names(bin)[names(bin) %in% var_index]
  out[1,] <- bin[1,names(bin) %in% var_index]
  out[2,] <- bin[2,names(bin) %in% var_index]
  out[3,] <- bin[3,names(bin) %in% var_index]
  
  out <- out[,var_index]
  return(out)
}

coef_plot.m3 <- get_point_estimates(m3, m3coeffs_cl, c("ag_opinion_submitted1",
                                                       "z.subject_matter.count",
                                                       "z.ms_pos_auto0","z.ms_pos_auto1",
                                                       "com_pos_auto.factorNo restrictions","com_pos_auto.factorRestrictions",
                                                       "panel_size","number_of_observations","legact_primary","legact_secondary"))


# plot logistic regression coefficients for figure 3
pdf(file = "figures/coefficients_ag.pdf", height = 6, width = 10)
par(mar=c(5,18,3,3))
plot("", axes = FALSE,
     xlim = c(-1.75,1.75),
     ylim = c(0.75,10.25),
     xlab = "Coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Post-Nice Treaty: Negative reference (Model 3)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-1.75,1.75,.25), cex.axis = 1.2)
axis(2,at = seq(10,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = c("AG opinion requested",
               "Number of subject matters",
               "MS favouring no restrictions",
               "MS favouring restrictions",
               "Commission position: No restrictions",
               "Commission position: Restrictions",
               "Panel size",
               "Number of observations",
               "Issue concerns primary EU law",
               "Issue concerns secondary EU law"))
abline(h = seq(10,1,-1), lty = 2, lwd = .5, col = "grey")
# m3
points(coef_plot.m3[1,], seq(10,1,-1), pch = 19)
segments(t(coef_plot.m3[2,]), seq(10,1,-1),
         t(coef_plot.m3[3,]), seq(10,1,-1), lwd = 2)
abline(v=0, lty = 2)
dev.off()



# multi-level logistic regression for table 1 -----------------------------

options(mc.cores = parallel::detectCores())

# multi-level models
options(mc.cores = parallel::detectCores())

# estimate model 4
stan.fit1 <- stan_glmer(mentioned_negative_dummy ~ observation_referring_ms + pos_auto.factor + oral_observation + panel_size + 
                          actor.collapsed + z.gdp + 
                          legact_primary + legact_secondary + z.subject_matter.count +
                          year_proceeding.factor + 
                          administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
                          labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport +
                          (1 | issue_id) + (1 | iuropa_proceeding_id),
                        data = mentions_submission_level, family = binomial(link = "logit"),
                        prior = normal(0,2.5), prior_intercept = normal(0,10),
                        chains=4, seed=123, iter=21000, warmup=1000)

save(stan.fit1, file = "objects/stan1_object.RData")
#load("objects/stan1_object.RData")
summary(stan.fit1)
ranefs.1 <- ranef(stan.fit1)

# get full results for fixed effect coefficients
stan1_id <- which(names(fixef(stan.fit1)) %in% c("observation_referring_ms",
                                                 "pos_auto.factorNo restrictions",
                                                 "pos_auto.factorRestrictions",
                                                 "oral_observation",
                                                 "z.subject_matter.count",
                                                 "panel_size",
                                                 "legact_primary",
                                                 "legact_secondary",
                                                 "z.gdp"))

BETA.1 <- as.matrix(stan.fit1)[,stan1_id]
mean.BETA.1 <- apply(BETA.1,2,mean)
sd.BETA.1 <- apply(BETA.1,2,sd)
low.BETA.1 <- apply(BETA.1,2,quantile,c(0.025,0.975))[1,]
high.BETA.1 <- apply(BETA.1,2,quantile,c(0.025,0.975))[2,]

# estimates for model 4 in table 1
round(mean.BETA.1, digits = 2)
round(low.BETA.1, digits = 2)
round(high.BETA.1, digits = 2)
var(unlist(ranef(stan.fit1)))


# estimate model 5
stan.fit2 <- stan_glmer(mentioned_negative_dummy ~ observation_referring_ms + pos_auto.factor + oral_observation + panel_size + 
                          actor.collapsed + z.gdp + 
                          legact_primary + legact_secondary + z.subject_matter.count +
                          year_proceeding.factor + 
                          administrative + agriculture + criminal + customs + environment + healthcare + intellectual_property +
                          labour + marketing + migration + procedure + procurement + no_area_of_law + social + tax + transport +
                          (1 | issue_id) + (1 | iuropa_proceeding_id),
                        data = subset_submissions_disagree, family = binomial(link = "logit"),
                        prior = normal(0,2.5), prior_intercept = normal(0,10),
                        chains=4, seed=123, iter=21000, warmup=1000)
summary(stan.fit2)
save(stan.fit2, file = "objects/stan2_object.RData")
#load("objects/stan2_object.RData")

# get full results for fixed effect coefficients
stan2_id <- which(names(fixef(stan.fit2)) %in% c("observation_referring_ms",
                                                 "pos_auto.factorNo restrictions",
                                                 "pos_auto.factorRestrictions",
                                                 "oral_observation",
                                                 "z.subject_matter.count",
                                                 "panel_size",
                                                 "legact_primary",
                                                 "legact_secondary",
                                                 "z.gdp"))
BETA.2 <- as.matrix(stan.fit2)[,stan2_id]
mean.BETA.2 <- apply(BETA.2,2,mean)
sd.BETA.2 <- apply(BETA.2,2,sd)
low.BETA.2 <- apply(BETA.2,2,quantile,c(0.025,0.975))[1,]
high.BETA.2 <- apply(BETA.2,2,quantile,c(0.025,0.975))[2,]

# estimates for model 5 in table 1
round(mean.BETA.2, digits = 2)
round(low.BETA.2, digits = 2)
round(high.BETA.2, digits = 2)
var(unlist(ranef(stan.fit2)))



# path analysis for figure 4 ----------------------------------------------

path3 <- '
  cjeu_pos_auto.factor ~ pos_auto.numerical + ag_pos_auto.numerical + com_pos_auto.numerical
  mentioned_negative_dummy ~ cjeu_pos_auto.factor + observation_referring_ms + pos_auto.numerical + 
    oral_observation + z.subject_matter.count + panel_size +
    legact_primary + legact_secondary + z.gdp + largest_ms +
    year_proceeding.factor
'
path3.fit <- sem(path3, data = mentions_submission_level)
summary(path3.fit)

# fit statistics for table 6 in the manuscript's supplementary materials
round(fitMeasures(path3.fit, c("chisq", "tli", "cfi", "rmsea", "pvalue")), 2)

GAMMA3_coefs <- parameterEstimates(path3.fit, ci = TRUE, level = 0.95)[1:3,c(3:4,8:9)]
BETA3_coefs <- parameterEstimates(path3.fit, ci = TRUE, level = 0.95)[4:13,c(3:4,8:9)]

# plot coefficients for GAMMA and BETA matrices for figure 4
pdf(file = "figures/submission_path.pdf", height = 8, width = 10)
par(mar=c(5,16,3,3), mfrow = c(2,1))
# GAMMA
plot("", axes = FALSE,
     xlim = c(-0.25,1.25),
     ylim = c(0.75,3.25),
     xlab = "Coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Path analysis: CJEU position (Model 6)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-0.25,1.25,.25), cex.axis = 1.2)
axis(2,at = seq(3,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = c("Member State position",
               "AG position",
               "Commission position"
     ))
abline(h = seq(3,1,-1), lty = 2, lwd = .5, col = "grey")
# path 3
points(GAMMA3_coefs[,2], seq(3,1,-1), pch = 19)
segments(GAMMA3_coefs[,3], seq(3,1,-1),
         GAMMA3_coefs[,4], seq(3,1,-1), lwd = 2)
abline(v=0, lty = 2)
# BETA
plot("", axes = FALSE,
     xlim = c(-0.05,0.11),
     ylim = c(0.75,10.25),
     xlab = "Coefficient estimates with 95% confidence intervals",
     ylab = "", pch = 19, cex = 1.5,
     main = "Path analysis: Negative reference (Model 6)", cex.lab = 1.2,
     frame = TRUE)
axis(1,at = seq(-0.05,0.11,.05), cex.axis = 1.2)
axis(2,at = seq(10,1,-1), las = 1, tick = TRUE, cex.axis = 1.2,
     label = c("CJEU position",
               "Case originated in MS",
               "Member State position",
               "Oral hearing",
               "Number of subject matters",
               "Panel size",
               "Issue concerns primary EU law",
               "Issue concerns secondary EU law",
               "Member State GDP per capita",
               "Largest Member States"))
abline(h = seq(10,1,-1), lty = 2, lwd = .5, col = "grey")
# path 3
points(BETA3_coefs[,2], c(seq(10,1,-1)), pch = 19)
segments(BETA3_coefs[,3], c(seq(10,1,-1)),
         BETA3_coefs[,4], c(seq(10,1,-1)), lwd = 2)
abline(v=0, lty = 2)
dev.off()
